Última actualización: 15 mayo 2020
La calibración del modelo se hace de forma conjunta en términos de cantidad y calidad de agua utilizando una función objetivo que indica el grado de similitud entre los valores simulados y observados de caudal diario (m3/s), la concentración diaria de nitrato (NO3) y la concentración diaria de fósforo total (PT). Para establecer las concentraciones diarias de PT y NO3 se utilizaron las relaciones exponenciales que se presentan en la sección Sistema de Monitoreo:muestreo de fósforo_y_nitrato
Para ello se utilizaron los siguientes paquetes de R: hydroPSO y SWATplusR.
El código utilizado para la calibración se muestra en el Anexo.
La función objetivo utiliza la eficiencia de Kling-Gupta (KGE). Se aplica a: caudal, fósforo y nitrato. Luego se añade una eficiencia relativa al sesgo de fósforo y nitrato. Finalmente se aplican pesos para dar mas importancia a el caudal así como a la parte baja de la cuenca. Esto se hace por dos razones:
La funcion objetivo (FO) se define por:
Fobj = 0.4Fobj|x = 23 + 0.6Fobj|x = 41
Donde 0.4 y 0.6 son los pesos que se aplican a las funciones objetivo (Fobj) en cada punto de control (x). Para Paso de los Troncos x = 23 y Paso Pache x = 41.
Fobj|x se define de la siguiente manera:
Fobj|x = 0.5KGE(Flow)|x + 0.1KGE(PT)|x + 0.1KGE(NO3)|x + 0.15PBE(PT)|x + 0.15PBE(NO3)|x
Donde KGE es la eficiencia de Kling-Gupta, PBE es una eficiencia relativa al sesgo y los coeficientes 0.5, 0.1 y 0.15 son los pesos utilizados.
$$KGE = 1- \sqrt{(r-1)^2 + (\frac{\sigma_{sim}}{\sigma_{obs}}-1)^2 + (\frac{\mu_{sim}}{\mu_{obs}}-1)^2}$$
r es la correlacion lineal entre las observaciones y las simulaciones. σobs es la desviación estándar en las observaciones. σsim la desviación estándar en las simulaciones. μobs es la media de las observaciones, μsim el promedio de las simulaciones. Esto fue calculado utilizando la función KGE del paquete hydroGOF.
Luego, PBE se define como uno menos el valor absoluto de la suma de los residuos dividida entre la suma de las observaciones.
$$PBE = 1 - abs\left(\frac{\sum(sim-obs)}{\sum(obs)}\right)$$
Se consideraron 22 parámetros en cada region. La cuanca alta está definida por las subcuencas “1,8,10,14,15,16,17,18,19,20,22” y la cuanca baja por “2,3,4,5,6,7,9,11,12,13,21,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41”. El rango de variación de cada parámetro así como el tipo de cambio se muestra en la siguiente tabla:
par_all <- readRDS("~/R_projects/SWAT_calibration_ETmodis/resultados/optim20200401_20/par_all.RDS")
best_sim <- readRDS("~/R_projects/SWAT_calibration_ETmodis/resultados/optim20200401_20/best_sim.RDS")
par_upr = apply(par_all,2,max)
par_lwr = apply(par_all,2,min)
par_upr = par_upr[best_sim$parameter$definition$par_name]
par_lwr = par_lwr[best_sim$parameter$definition$par_name]
par_optim = data.frame(parametro=best_sim$parameter$definition$par_name,
archivo=best_sim$parameter$definition$file_name,
cambio=best_sim$parameter$definition$change,
L_inf = as.numeric(round(par_lwr,2)),
L_sup = as.numeric(round(par_upr,2)),
Best_par = as.numeric(round(best_sim$parameter$values,2)))
kable(par_optim, caption ="Tipo de cambio, rangos de los parámetros en la optimización y mejor parámetro") %>%
kable_styling(c("striped", "bordered"))| parametro | archivo | cambio | L_inf | L_sup | Best_par |
|---|---|---|---|---|---|
| r1_CN2 | mgt | pctchg | -0.40 | 0.00 | -38.42 |
| r1_ESCO | hru | absval | 0.10 | 0.40 | 0.25 |
| r1_SOL_AWC | sol | pctchg | 0.00 | 0.20 | 16.44 |
| r1_ALPHA_BF | gw | absval | 0.40 | 0.80 | 0.52 |
| r1_GWQMN | gw | absval | 800.00 | 6000.00 | 3487.58 |
| r1_GW_DELAY | gw | absval | 20.00 | 35.00 | 22.65 |
| r1_REVAPMN | gw | absval | 325.00 | 375.00 | 347.13 |
| r1_GW_REVAP | gw | absval | 0.02 | 0.10 | 0.07 |
| r1_OV_N | hru | absval | 0.40 | 0.60 | 0.53 |
| r1_SLSUBBSN | hru | pctchg | -0.20 | 0.05 | -11.31 |
| r2_CN2 | mgt | pctchg | -0.20 | 0.20 | 8.48 |
| r2_ESCO | hru | absval | 0.35 | 0.45 | 0.40 |
| r2_SOL_AWC | sol | pctchg | -0.15 | 0.00 | -4.76 |
| r2_ALPHA_BF | gw | absval | 0.40 | 0.80 | 0.58 |
| r2_GWQMN | gw | absval | 800.00 | 1500.00 | 890.45 |
| r2_GW_DELAY | gw | absval | 25.00 | 45.00 | 41.10 |
| r2_REVAPMN | gw | absval | 610.00 | 710.00 | 695.42 |
| r2_GW_REVAP | gw | absval | 0.03 | 0.06 | 0.05 |
| r2_OV_N | hru | absval | 0.40 | 0.60 | 0.45 |
| r2_SLSUBBSN | hru | pctchg | 0.20 | 0.40 | 29.61 |
| r1_HRU_SLP | hru | pctchg | -0.10 | 0.10 | -4.08 |
| r1_USLE_K | sol | pctchg | -0.15 | 0.15 | -0.42 |
| r1_USLE_P | mgt | pctchg | -0.15 | 0.15 | -10.36 |
| r1_P_UPDIS | bsn | absval | 20.00 | 60.00 | 43.82 |
| r1_FRT_SURFACE | mgt | absval | 0.10 | 0.60 | 0.37 |
| r1_FRT_KG | mgt | pctchg | -0.45 | -0.15 | -34.07 |
| r1_NPERCO | bsn | absval | 0.50 | 0.80 | 0.68 |
| r1_ERORGN | hru | absval | 3.00 | 4.00 | 3.33 |
| r2_HRU_SLP | hru | pctchg | -0.10 | 0.10 | -8.21 |
| r2_USLE_K | sol | pctchg | 0.00 | 0.20 | 0.87 |
| r2_USLE_P | mgt | pctchg | -0.15 | 0.15 | -4.84 |
| r2_P_UPDIS | bsn | absval | 40.00 | 90.00 | 78.00 |
| r2_FRT_SURFACE | mgt | absval | 0.10 | 0.60 | 0.26 |
| r2_FRT_KG | mgt | pctchg | -0.40 | 0.10 | -30.88 |
| r2_NPERCO | bsn | absval | 0.70 | 1.00 | 0.98 |
| r2_ERORGN | hru | absval | 2.00 | 4.00 | 2.13 |
| r1_PPERCO | bsn | absval | 10.00 | 17.50 | 14.63 |
| r2_PPERCO | bsn | absval | 10.00 | 17.50 | 15.59 |
| r1_CDN | bsn | absval | 0.00 | 3.00 | 0.79 |
| r2_CDN | bsn | absval | 0.00 | 3.00 | 2.61 |
| r1_CMN | bsn | absval | 0.00 | 0.03 | 0.02 |
| r2_CMN | bsn | absval | 0.00 | 0.03 | 0.02 |
| r1_PSP | bsn | absval | 0.01 | 0.80 | 0.43 |
| r2_PSP | bsn | absval | 0.01 | 0.80 | 0.15 |
Se muestra el KGE a escala diaria de las mejores simulaciones.
best_skill <- readRDS("~/R_projects/SWAT_calibration_ETmodis/resultados/optim20200401_20/best_skill.RDS")
library(plotly)
fig <- plot_ly(y = ~best_skill[,"KGE41"], type = "box",name="Paso Pache \n(calibración)" )
fig <- fig %>% add_trace(y = ~best_skill[,"KGE36"],name="Fray Marcos \n(validación)")
fig <- fig %>% add_trace(y = ~best_skill[,"KGE23"],name="Paso de los Troncos \n(calibración)")
fig <- fig %>% layout(yaxis = list(title="KGE [-]"), showlegend = FALSE)
figLas simulaciones en Paso de los Troncos sobrestiman el caudal.
fig <- plot_ly(y = ~best_skill[,"pbias41"], type = "box",name="Paso Pache \n(calibración)" )
fig <- fig %>% add_trace(y = ~best_skill[,"pbias36"],name="Fray Marcos \n(validación)")
fig <- fig %>% add_trace(y = ~best_skill[,"pbias23"],name="Paso de los Troncos \n(calibración)")
fig <- fig %>% layout(yaxis = list(title="BIAS [%]"), showlegend = FALSE)
figEn paso de los troncos se sobrestima el caudal especialmente después de 2013.
library(zoo)
obs.data <- readRDS("~/R_projects/SWAT_calibration_ETmodis/data/obs.data.RDS")
obs.data = data.frame(obs.data)
obs.data$q23lag1 = c(obs.data$q23[2:nrow(obs.data)],NA)
obs.data$q36lag1 = c(obs.data$q36[2:nrow(obs.data)],NA)
obs.data$q41lag2 = c(obs.data$q41[3:nrow(obs.data)],NA,NA)
Qobs = obs.data$q23lag1
Qsim = best_sim$simulation$FLOW_23
Rdate = best_sim$simulation$date
fig <- plot_ly(x = ~Rdate, y = ~Qobs, type="scatter",mode = "lines",
name="Qobs", line = list(color = 'black', width = 2))
fig <- fig %>% add_lines(x = ~Rdate, y = ~Qsim, name="Qsim",
type = 'scatter',
mode = 'line',
line = list(color = 'orange', width = 2))
fig <- fig %>% layout(legend = list(x = 0.1, y = 0.9),yaxis = list(title="Flow (m3/s)"),xaxis = list(title="") )
figQobs = obs.data$q36lag1
Qsim = best_sim$simulation$FLOW_36
Rdate = best_sim$simulation$date
fig <- plot_ly(x = ~Rdate, y = ~Qobs, type="scatter",mode = "lines",
name="Qobs", line = list(color = 'black', width = 2))
fig <- fig %>% add_lines(x = ~Rdate, y = ~Qsim, name="Qsim",
type = 'scatter',
mode = 'line',
line = list(color = 'orange', width = 2))
fig <- fig %>% layout(legend = list(x = 0.1, y = 0.9),yaxis = list(title="Flow (m3/s)"),xaxis = list(title="") )
figQobs = obs.data$q41lag2
Qsim = best_sim$simulation$FLOW_41
Rdate = best_sim$simulation$date
fig <- plot_ly(x = ~Rdate, y = ~Qobs, type="scatter",mode = "lines",
name="Qobs", line = list(color = 'black', width = 2))
fig <- fig %>% add_lines(x = ~Rdate, y = ~Qsim, name="Qsim",
type = 'scatter',
mode = 'line',
line = list(color = 'orange', width = 2))
fig <- fig %>% layout(legend = list(x = 0.1, y = 0.9),yaxis = list(title="Flow (m3/s)"),xaxis = list(title="") )
figLos resultados de fósforo total en Paso Pache son bastante buenos, no obstante algo no va bien con las simulaciones en Paso de los Troncos. Posiblemente eso sea un problema heredado de una mala calibración en caudal.
# datos observados
area_41 = 4896
area_36 = 2744
area_23 = 687
qlt_41 <- readRDS("~/R_projects/SWAT_calibration_ETmodis/data/rch41_calidad.RDS")
qlt_23 <- readRDS("~/R_projects/SWAT_calibration_ETmodis/data/rch23_calidad.RDS")
qlt_41$qobs = 1000*qlt_41$qobs/area_41
qlt_23$qobs = 1000*qlt_23$qobs/area_23
qlt_41$NO3_mgxs = qlt_41$NO3_mgxs/area_41
qlt_23$NO3_mgxs = qlt_23$NO3_mgxs/area_23
qlt_41$PT_mgxs = qlt_41$PT_mgxs/area_41
qlt_23$PT_mgxs = qlt_23$PT_mgxs/area_23
fit41_PT = lm(log(qlt_41$PT_mgxs)~log(qlt_41$qobs))
fit41_NO3 = lm(log(qlt_41$NO3_mgxs)~log(qlt_41$qobs))
fit23_PT = lm(log(qlt_23$PT_mgxs)~log(qlt_23$qobs))
fit23_NO3 = lm(log(qlt_23$NO3_mgxs)~log(qlt_23$qobs))
fit41_PT = round(c(exp(fit41_PT$coefficients[1]), fit41_PT$coefficients[2]),2)
fit41_NO3 = round(c(exp(fit41_NO3$coefficients[1]), fit41_NO3$coefficients[2]),2)
fit23_PT = round(c(exp(fit23_PT$coefficients[1]), fit23_PT$coefficients[2]),2)
fit23_NO3 = round(c(exp(fit23_NO3$coefficients[1]), fit23_NO3$coefficients[2]),2)
qlt_41$fit_PT = fit41_PT[1]*qlt_41$qobs^fit41_PT[2]
qlt_41$fit_NO3 = fit41_NO3[1]*qlt_41$qobs^fit41_NO3[2]
qlt_23$fit_PT = fit23_PT[1]*qlt_23$qobs^fit23_PT[2]
qlt_23$fit_NO3 = fit23_NO3[1]*qlt_23$qobs^fit23_NO3[2]
# datos simulados
qsim = best_sim$simulation
qlt.sim_41 = data.frame(flow_lps_km2=1000*qsim$FLOW_41/4896,
PT_mgs_km2 = 1000000*qsim$PT_41/(24*60*60*4896),
NO3_mgs_km2 = 1000000*qsim$NO3_41/(24*60*60*4896),
qobs_lps_km2 = 1000*obs.data$q41lag2/4896,
PTobs_mgs_km2 = 0.18*(1000*obs.data$q41lag2/4896)^1.02,
NO3obs_mgs_km2 = 0.11*(1000*obs.data$q41lag2/4896)^1.41)
qlt.sim_41[qlt.sim_41<=0] = NA
qlt.sim_41 = na.omit(qlt.sim_41)
qlt.sim_23 = data.frame(flow_lps_km2=1000*qsim$FLOW_23/687,
PT_mgs_km2 = 1000000*qsim$PT_23/(24*60*60*687),
NO3_mgs_km2 = 1000000*qsim$NO3_23/(24*60*60*687),
qobs_lps_km2 = 1000*obs.data$q23lag1/687,
PTobs_mgs_km2 = 0.1*(1000*obs.data$q23lag1/687)^0.93,
NO3obs_mgs_km2 = 0.08*(1000*obs.data$q23lag1/687)^1.41)
qlt.sim_23[qlt.sim_23<=0] = NA
qlt.sim_23 = na.omit(qlt.sim_23)
fit = lm(log(PT_mgs_km2)~log(flow_lps_km2), data=qlt.sim_41)
a = exp(fit$coefficients[1])
b = fit$coefficients[2]
q41fit = seq(min(qlt.sim_41$flow_lps_km2), max(qlt.sim_41$flow_lps_km2), length.out = 30)
PT41fit = a*q41fit^b
PT = plot_ly(data = qlt.sim_41, x=~flow_lps_km2,y=~PT_mgs_km2, type='scatter',
mode='markers', name="Simulado",
marker=list(color='gray', size=4))
PT = PT %>% add_markers(data=qlt_41,x = ~qobs, y = ~PT_mgxs,
type = 'scatter',
mode = 'markers',
marker = list(color = 'red', size=8),
name="Observado")
PT <- PT %>% add_lines(data=qlt_41,x = ~qobs, y = ~fit_PT,
type = 'scatter',
mode = 'line',
marker = list(color = 'red', size=0),
line = list(color = 'red', width = 2),
name="Regresión observados")
PT <- PT %>% add_lines(x= ~q41fit, y = ~PT41fit,
type = 'scatter',
mode = 'line',
marker = list(color = 'black', size=0),
line = list(color = 'black', width = 2),
name="Regresión simulados")
PT <- layout(PT, xaxis = list(title="Caudal (lps/km2)",range=c(-1,2),type = "log"),
yaxis = list(title="PT (mg/s/km2)", range=c(-2,1),type = "log"),
legend = list(x = 0.01, y = 1)
# showlegend = FALSE,
)
PTfit = lm(log(PT_mgs_km2)~log(flow_lps_km2), data=qlt.sim_23)
a = exp(fit$coefficients[1])
b = fit$coefficients[2]
q23fit = seq(min(qlt.sim_23$flow_lps_km2), max(qlt.sim_23$flow_lps_km2), length.out = 30)
PT23fit = a*q23fit^b
PT = plot_ly(data = qlt.sim_23, x=~flow_lps_km2,y=~PT_mgs_km2, type='scatter',
mode='markers', name="Simulado",
marker=list(color='gray', size=4))
PT = PT %>% add_markers(data=qlt_23,x = ~qobs, y = ~PT_mgxs,
type = 'scatter',
mode = 'markers',
marker = list(color = 'red', size=8),
name="Observado")
PT <- PT %>% add_lines(data=qlt_23,x = ~qobs, y = ~fit_PT,
type = 'scatter',
mode = 'line',
marker = list(color = 'red', size=0),
line = list(color = 'red', width = 2),
name="Regresión observados")
PT <- PT %>% add_lines(x= ~q23fit, y = ~PT23fit,
type = 'scatter',
mode = 'line',
marker = list(color = 'black', size=1),
line = list(color = 'black', width = 2),
name="Regresión simulados")
PT <- layout(PT, xaxis = list(title="Caudal (lps/km2)",range=c(-1,2),type = "log"),
yaxis = list(title="PT (mg/s/km2)", range=c(-2,1),type = "log"),
legend = list(x = 0.01, y = 1)
# showlegend = FALSE,
)
PTfit = lm(log(NO3_mgs_km2)~log(flow_lps_km2), data=qlt.sim_41)
a = exp(fit$coefficients[1])
b = fit$coefficients[2]
q41fit = seq(min(qlt.sim_41$flow_lps_km2), max(qlt.sim_41$flow_lps_km2), length.out = 30)
NO3_41fit = a*q41fit^b
PT = plot_ly(data = qlt.sim_41, x=~flow_lps_km2,y=~NO3_mgs_km2, type='scatter',
mode='markers', name="Simulado",
marker=list(color='gray', size=4))
PT = PT %>% add_markers(data=qlt_41,x = ~qobs, y = ~NO3_mgxs,
type = 'scatter',
mode = 'markers',
marker = list(color = 'red', size=8),
name="Observado")
PT <- PT %>% add_lines(data=qlt_41,x = ~qobs, y = ~fit_NO3,
type = 'scatter',
mode = 'line',
marker = list(color = 'red', size=0),
line = list(color = 'red', width = 2),
name="Regresión observados")
PT <- PT %>% add_lines(x= ~q41fit, y = ~NO3_41fit,
type = 'scatter',
mode = 'line',
marker = list(color = 'black', size=1),
line = list(color = 'black', width = 2),
name="Regresión simulados")
PT <- layout(PT, xaxis = list(title="Caudal (lps/km2)",range=c(-1,2),type = "log"),
yaxis = list(title="NO3 (mg/s/km2)", range=c(-2,1),type = "log"),
legend = list(x = 0.7, y = 0.01)
# showlegend = FALSE,
)
PTfit = lm(log(NO3_mgs_km2)~log(flow_lps_km2), data=qlt.sim_23)
a = exp(fit$coefficients[1])
b = fit$coefficients[2]
q23fit = seq(min(qlt.sim_23$flow_lps_km2), max(qlt.sim_23$flow_lps_km2), length.out = 30)
NO3_23fit = a*q23fit^b
PT = plot_ly(data = qlt.sim_23, x=~flow_lps_km2,y=~NO3_mgs_km2, type='scatter',
mode='markers', name="Simulado",
marker=list(color='gray', size=4))
PT = PT %>% add_markers(data=qlt_23,x = ~qobs, y = ~NO3_mgxs,
type = 'scatter',
mode = 'markers',
marker = list(color = 'red', size=8),
name="Observado")
PT <- PT %>% add_lines(data=qlt_23,x = ~qobs, y = ~fit_NO3,
type = 'scatter',
mode = 'line',
marker = list(color = 'red', size=0),
line = list(color = 'red', width = 2),
name="Regresión observados")
PT <- PT %>% add_lines(x= ~q23fit, y = ~NO3_23fit,
type = 'scatter',
mode = 'line',
marker = list(color = 'black', size=1),
line = list(color = 'black', width = 2),
name="Regresión simulados")
PT <- layout(PT, xaxis = list(title="Caudal (lps/km2)",range=c(-1,2),type = "log"),
yaxis = list(title="NO3 (mg/s/km2)", range=c(-2,1),type = "log"),
legend = list(x = 0.6, y = 0.01)
# showlegend = FALSE,
)
PTbest_par <- readRDS("~/R_projects/SWAT_calibration_ETmodis/resultados/optim20200401_20/best_par.RDS")
par_name = names(par_lwr)[order(names(par_lwr))]
par(mfrow=c(12,2), mar=c(4,4,1,1))
for(i in 1:22){
p1 = par_name[i]
p2 = par_name[i+22]
hist(best_par[,p1], main=paste(p1),
breaks = seq(par_lwr[p1], par_upr[p1], length.out = 30))
hist(best_par[,p2], main=paste(p2),
breaks = seq(par_lwr[p2], par_upr[p2], length.out = 30))
}rm(list = ls())
# Librerias
library(SWATplusR)
library(stringr)
library(hydroGOF)
library(hydromad)
library(tibble)
library(hydroPSO)
# Caudales observados
obs.data <- readRDS("./data/obs.data.RDS")
obs.data = as.data.frame(obs.data)
obs.data$q23lag1 = c(obs.data$q23[2:nrow(obs.data)],NA)
obs.data$q36lag1 = c(obs.data$q36[2:nrow(obs.data)],NA)
obs.data$q41lag2 = c(obs.data$q41[3:nrow(obs.data)],NA,NA)
obs.data$PT_23_lag1 = c(obs.data$PT_23[2:nrow(obs.data)],NA)
obs.data$NO3_23_lag1 = c(obs.data$NO3_23[2:nrow(obs.data)],NA)
obs.data$PT_41_lag2 = c(lag(obs.data$PT_41)[3:nrow(obs.data)],NA,NA)
obs.data$NO3_41_lag2 = c(lag(obs.data$NO3_41)[3:nrow(obs.data)],NA,NA)
# SWAT project mensual
swat_project = "./SWAT_SL_20200401/"
# Funcion de optimización
swat_optim <- function(mypar, obs.data, project_dir, start_sim, end_sim,runrun_path, head_log) {
path_run = paste0(runrun_path,"run_",Sys.getpid())
par_set = tibble("r1_CN2::CN2.mgt|change = pctchg |
sub == c(1,8,10,14,15,16,17,18,19,20,22)" =
mypar[1]*100,
"r1_ESCO::ESCO.hru|change = absval |
sub == c(1,8,10,14,15,16,17,18,19,20,22)" =
mypar[2],
"r1_SOL_AWC::SOL_AWC.sol|change = pctchg |
sub == c(1,8,10,14,15,16,17,18,19,20,22)" =
mypar[3]*100,
"r1_ALPHA_BF::ALPHA_BF.gw|change = absval |
sub == c(1,8,10,14,15,16,17,18,19,20,22)" =
mypar[4],
"r1_GWQMN::GWQMN.gw|change = absval |
sub == c(1,8,10,14,15,16,17,18,19,20,22)" =
mypar[5],
"r1_GW_DELAY::GW_DELAY.gw|change = absval |
sub == c(1,8,10,14,15,16,17,18,19,20,22)" =
mypar[6],
"r1_REVAPMN::REVAPMN.gw|change = absval |
sub == c(1,8,10,14,15,16,17,18,19,20,22)" =
mypar[7],
"r1_GW_REVAP::GW_REVAP.gw|change = absval |
sub == c(1, 16, 18, 22,23)" =
mypar[8],
"r1_OV_N::OV_N.hru|change = absval |
sub == c(1,8,10,14,15,16,17,18,19,20,22)" =
mypar[9],
"r1_SLSUBBSN::SLSUBBSN.hru|change = pctchg |
sub == c(1,8,10,14,15,16,17,18,19,20,22)" =
mypar[10]*100,
"r2_CN2::CN2.mgt|change = pctchg |
sub == c(2,3,4,5,6,7,9,11,12,13,21,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41)" =
mypar[11]*100,
"r2_ESCO::ESCO.hru|change = absval |
sub == c(2,3,4,5,6,7,9,11,12,13,21,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41)" =
mypar[12],
"r2_SOL_AWC::SOL_AWC.sol|change = pctchg |
sub == c(2,3,4,5,6,7,9,11,12,13,21,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41)" =
mypar[13]*100,
"r2_ALPHA_BF::ALPHA_BF.gw|change = absval |
sub == c(2,3,4,5,6,7,9,11,12,13,21,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41)" =
mypar[14],
"r2_GWQMN::GWQMN.gw|change = absval |
sub == c(2,3,4,5,6,7,9,11,12,13,21,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41)" =
mypar[15],
"r2_GW_DELAY::GW_DELAY.gw|change = absval |
sub == c(2,3,4,5,6,7,9,11,12,13,21,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41)" =
mypar[16],
"r2_REVAPMN::REVAPMN.gw|change = absval |
sub == c(2,3,4,5,6,7,9,11,12,13,21,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41)" =
mypar[17],
"r2_GW_REVAP::GW_REVAP.gw|change = absval |
sub == c(2,3,4,5,6,7,9,11,12,13,21,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41)" =
mypar[18],
"r2_OV_N::OV_N.hru|change = absval |
sub == c(2,3,4,5,6,7,9,11,12,13,21,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41)" =
mypar[19],
"r2_SLSUBBSN::SLSUBBSN.hru|change = pctchg |
sub == c(2,3,4,5,6,7,9,11,12,13,21,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41)" =
mypar[20]*100,
"r1_HRU_SLP::HRU_SLP.hru|change= pctchg |
sub == c(1,8,10,14,15,16,17,18,19,20,22)"= mypar[21]*100,
"r1_USLE_K::USLE_K.sol|change= pctchg |
sub == c(1,8,10,14,15,16,17,18,19,20,22)" = mypar[22]*100,
"r1_USLE_P::USLE_P.mgt|change= pctchg |
sub == c(1,8,10,14,15,16,17,18,19,20,22)"=mypar[23]*100,
"r1_P_UPDIS::P_UPDIS.bsn|change= absval |
sub == c(1,8,10,14,15,16,17,18,19,20,22)"= mypar[24],
"r1_FRT_SURFACE::FRT_SURFACE.mgt|change=absval |
sub == c(1,8,10,14,15,16,17,18,19,20,22)" = mypar[25],
"r1_FRT_KG::FRT_KG.mgt|change=pctchg |
sub == c(1,8,10,14,15,16,17,18,19,20,22)" = mypar[26]*100,
"r1_NPERCO::NPERCO.bsn|change= absval |
sub == c(1,8,10,14,15,16,17,18,19,20,22)"=mypar[27],
"r1_ERORGN::ERORGN.hru|change= absval |
sub == c(1,8,10,14,15,16,17,18,19,20,22)"= mypar[28],
"r2_HRU_SLP::HRU_SLP.hru|change= pctchg |
sub == c(2,3,4,5,6,7,9,11,12,13,21,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41)"=
mypar[29]*100,
"r2_USLE_K::USLE_K.sol|change= pctchg |
sub == c(2,3,4,5,6,7,9,11,12,13,21,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41)" =
mypar[30]*100,
"r2_USLE_P::USLE_P.mgt|change= pctchg |
sub == c(2,3,4,5,6,7,9,11,12,13,21,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41)"=
mypar[31]*100,
"r2_P_UPDIS::P_UPDIS.bsn|change= absval |
sub == c(2,3,4,5,6,7,9,11,12,13,21,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41)"=
mypar[32],
"r2_FRT_SURFACE::FRT_SURFACE.mgt|change=absval |
sub == c(2,3,4,5,6,7,9,11,12,13,21,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41)"=
mypar[33],
"r2_FRT_KG::FRT_KG.mgt|change=pctchg |
sub == c(2,3,4,5,6,7,9,11,12,13,21,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41)"=
mypar[34]*100,
"r2_NPERCO::NPERCO.bsn|change= absval |
sub == c(2,3,4,5,6,7,9,11,12,13,21,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41)"=
mypar[35],
"r2_ERORGN::ERORGN.hru|change= absval |
sub == c(2,3,4,5,6,7,9,11,12,13,21,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41)"=
mypar[36],
"r1_PPERCO::PPERCO.bsn|change= absval |
sub == c(1,8,10,14,15,16,17,18,19,20,22)"= mypar[37],
"r2_PPERCO::PPERCO.bsn|change= absval |
sub == c(2,3,4,5,6,7,9,11,12,13,21,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41)"=
mypar[38],
"r1_CDN::CDN.bsn|change= absval |
sub == c(1,8,10,14,15,16,17,18,19,20,22)"= mypar[39],
"r2_CDN::CDN.bsn|change= absval |
sub == c(2,3,4,5,6,7,9,11,12,13,21,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41)"=
mypar[40],
"r1_CMN::CMN.bsn|change= absval |
sub == c(1,8,10,14,15,16,17,18,19,20,22)"= mypar[41],
"r2_CMN::CMN.bsn|change= absval |
sub == c(2,3,4,5,6,7,9,11,12,13,21,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41)"=
mypar[42],
"r1_PSP::PSP.bsn|change= absval |
sub == c(1,8,10,14,15,16,17,18,19,20,22)"= mypar[43],
"r2_PSP::PSP.bsn|change= absval |
sub == c(2,3,4,5,6,7,9,11,12,13,21,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41)"=
mypar[44]
)
# run swat
qsim = run_swat2012(project_path = project_dir,
output = list(FLOW = define_output(file = "rch",
variable = "FLOW_OUT",
unit = c(23,36,41)),
PT = define_output(file = "rch",
variable = "TOT_Pkg",
unit = c(23,41)),
NO3 = define_output(file = "rch",
variable = "NO3_OUT",
unit = c(23,41))),
parameter = par_set,
start_date = as.Date("2009-01-01"),
end_date = as.Date("2015-12-31"),
years_skip = 1,
run_path=path_run,
quiet = TRUE,
output_interval="d",
add_parameter = FALSE)
# objective function
nse23 = NSE(as.vector(qsim$FLOW_23), as.vector(obs.data$q23lag1))
nse36 = NSE(as.vector(qsim$FLOW_36), as.vector(obs.data$q36lag1))
nse41 = NSE(as.vector(qsim$FLOW_41), as.vector(obs.data$q41lag2))
kge23 = KGE(as.vector(qsim$FLOW_23), as.vector(obs.data$q23lag1))
kge36 = KGE(as.vector(qsim$FLOW_36), as.vector(obs.data$q36lag1))
kge41 = KGE(as.vector(qsim$FLOW_41), as.vector(obs.data$q41lag2))
pbias23 = pbias(as.vector(qsim$FLOW_23), as.vector(obs.data$q23lag1))
pbias36 = pbias(as.vector(qsim$FLOW_36), as.vector(obs.data$q36lag1))
pbias41 = pbias(as.vector(qsim$FLOW_41), as.vector(obs.data$q41lag2))
# objective quality
qlt_41 = data.frame(flow_lps_km2=1000*qsim$FLOW_41/4896,
PT_mgs_km2 = 1000000*qsim$PT_41/(24*60*60*4896),
NO3_mgs_km2 = 1000000*qsim$NO3_41/(24*60*60*4896),
qobs_lps_km2 = 1000*obs.data$q41lag2/4896,
PTobs_mgs_km2 = 0.18*(1000*obs.data$q41lag2/4896)^1.02,
NO3obs_mgs_km2 = 0.11*(1000*obs.data$q41lag2/4896)^1.41)
qlt_41[qlt_41<=0] = NA
qlt_41 = na.omit(qlt_41)
qlt_23 = data.frame(flow_lps_km2=1000*qsim$FLOW_23/687,
PT_mgs_km2 = 1000000*qsim$PT_23/(24*60*60*687),
NO3_mgs_km2 = 1000000*qsim$NO3_23/(24*60*60*687),
qobs_lps_km2 = 1000*obs.data$q23lag1/687,
PTobs_mgs_km2 = 0.1*(1000*obs.data$q23lag1/687)^0.93,
NO3obs_mgs_km2 = 0.08*(1000*obs.data$q23lag1/687)^1.41)
qlt_23[qlt_23<=0] = NA
qlt_23 = na.omit(qlt_23)
fit = lm(log(PT_mgs_km2)~log(flow_lps_km2), data=qlt_41)
a = exp(fit$coefficients[1])
b = fit$coefficients[2]
qlt_41$PT_fit = a*(qlt_41$flow_lps_km2)^b
fit = lm(log(NO3_mgs_km2)~log(flow_lps_km2), data=qlt_41)
a = exp(fit$coefficients[1])
b = fit$coefficients[2]
qlt_41$NO3_fit = a*(qlt_41$flow_lps_km2)^b
fit = lm(log(PT_mgs_km2)~log(flow_lps_km2), data=qlt_23)
a = exp(fit$coefficients[1])
b = fit$coefficients[2]
qlt_23$PT_fit = a*(qlt_23$flow_lps_km2)^b
fit = lm(log(NO3_mgs_km2)~log(flow_lps_km2), data=qlt_23)
a = exp(fit$coefficients[1])
b = fit$coefficients[2]
qlt_23$NO3_fit = a*(qlt_23$flow_lps_km2)^b
kge_PT_41 = KGE(qlt_41$PT_fit, qlt_41$PTobs_mgs_km2)
kge_NO3_41 = KGE(qlt_41$NO3_fit, qlt_41$NO3obs_mgs_km2)
kge_PT_23 = KGE(qlt_23$PT_fit, qlt_23$PTobs_mgs_km2)
kge_NO3_23 = KGE(qlt_23$NO3_fit, qlt_23$NO3obs_mgs_km2)
pbias_PT_23 = pbias(qlt_23$PT_fit, qlt_23$PTobs_mgs_km2)
pbias_PT_41 = pbias(qlt_41$PT_fit, qlt_41$PTobs_mgs_km2)
pbias_NO3_23 = pbias(qlt_23$NO3_fit, qlt_23$NO3obs_mgs_km2)
pbias_NO3_41 = pbias(qlt_41$NO3_fit, qlt_41$NO3obs_mgs_km2)
FO23 = 0.5*kge23 + 0.1*kge_NO3_23 + 0.1*kge_PT_23 + 0.15*(1-abs(pbias_PT_23)/100) +
0.15*(1-abs(pbias_NO3_23)/100)
FO41 = 0.5*kge41 + 0.1*kge_NO3_41 + 0.1*kge_PT_41 + 0.15*(1-abs(pbias_PT_41)/100) +
0.15*(1-abs(pbias_NO3_41)/100)
FO = 0.6*FO41 + 0.4*FO23
# logfile
vout = data.frame(FO,kge23,kge36,kge41,
nse23,nse36,nse41,
pbias23,pbias36,pbias41,
kge_PT_23, kge_PT_41, kge_NO3_41, kge_NO3_23,
pbias_PT_23, pbias_PT_41, pbias_NO3_23, pbias_NO3_41,
t(mypar))
colnames(vout) = head_log
log.file = paste0(runrun_path,"log_",Sys.getpid(),".txt")
if(!any(list.files(runrun_path)==paste0("log_",Sys.getpid(),".txt"))){
write.table(t(head_log), log.file, col.names = FALSE, row.names = FALSE, quote = FALSE)
}
ff = read.table(log.file, col.names = head_log,colClasses = "character")
ff = rbind(ff,vout)
write.table(ff, log.file, col.names = FALSE, row.names = FALSE, quote = FALSE)
#output
return(FO)
}
r1_CN2 = c(-0.4, 0) ###
r2_CN2 = c(-0.2,0.2)
r1_ESCO = c(0.1,0.4)
r2_ESCO = c(0.35,0.45)
r1_SOL_AWC = c(0,0.2)
r2_SOL_AWC = c(-0.15,0)
r1_ALPHA_BF = c(0.4,0.8)
r2_ALPHA_BF = c(0.4,0.8)
r1_GWQMN = c(800,6000)
r2_GWQMN = c(800,1500)
r1_GW_DELAY = c(20,35)
r2_GW_DELAY = c(25,45)
r1_REVAPMN = c(325,375)
r2_REVAPMN = c(610,710)
r1_GW_REVAP = c(0.02, 0.1) ###
r2_GW_REVAP = c(0.03, 0.055)
r1_OV_N = c(0.4, 0.6)
r2_OV_N = c(0.4, 0.6)
r1_SLSUBBSN = c(-0.2, 0.05) ###
r2_SLSUBBSN = c(0.2, 0.4)
r1_HRU_SLP = c(-0.1,0.1)
r2_HRU_SLP = c(-0.1,0.1)
r1_USLE_K = c(-0.15,0.15)
r2_USLE_K = c(0,0.2)
r1_USLE_P = c(-0.15,0.15)
r2_USLE_P = c(-0.15,0.15)
r1_P_UPDIS = c(20, 60)
r2_P_UPDIS = c(40, 90)
r1_FRT_SURFACE = c(0.1,0.6)
r2_FRT_SURFACE = c(0.1,0.6)
r1_FRT_KG = c(-0.45,-0.15)
r2_FRT_KG = c(-0.4,0.1)
r1_NPERCO = c(0.5,0.8)
r2_NPERCO = c(0.7,1)
r1_ERORGN = c(3,4)
r2_ERORGN = c(2,4)
r1_PPERCO = c(10,17.5)
r2_PPERCO = c(10,17.5) #bsn absval
r1_CDN = c(0,3) #bsn absval
r2_CDN = c(0,3) #bsn absval
r1_CMN = c(0,0.03) #bsn absval
r2_CMN = c(0,0.03) #bsn absval
r1_PSP = c(0.01,0.8) #bsn absval
r2_PSP = c(0.01,0.8)
# CN2 ESCO SOL_AWC ALPHA_BF GWQMN GW_DELAY REVAPMN GW_REVAP OV_N SLSUBBSN
par_lwr = c(r1_CN2[1], r1_ESCO[1], r1_SOL_AWC[1],r1_ALPHA_BF[1], r1_GWQMN[1], r1_GW_DELAY[1],
r1_REVAPMN[1], r1_GW_REVAP[1], r1_OV_N[1], r1_SLSUBBSN[1],
r2_CN2[1], r2_ESCO[1], r2_SOL_AWC[1],r2_ALPHA_BF[1], r2_GWQMN[1], r2_GW_DELAY[1],
r2_REVAPMN[1], r2_GW_REVAP[1], r2_OV_N[1], r2_SLSUBBSN[1],
r1_HRU_SLP[1], r1_USLE_K[1], r1_USLE_P[1], r1_P_UPDIS[1], r1_FRT_SURFACE[1],
r1_FRT_KG[1], r1_NPERCO[1], r1_ERORGN[1],
r2_HRU_SLP[1], r2_USLE_K[1], r2_USLE_P[1], r2_P_UPDIS[1], r2_FRT_SURFACE[1],
r2_FRT_KG[1], r2_NPERCO[1], r2_ERORGN[1],
r1_PPERCO[1], r1_CDN[1], r1_CMN[1], r1_PSP[1],
r2_PPERCO[1], r2_CDN[1], r2_CMN[1], r2_PSP[1])
par_upr = c(r1_CN2[2], r1_ESCO[2], r1_SOL_AWC[2],r1_ALPHA_BF[2], r1_GWQMN[2], r1_GW_DELAY[2],
r1_REVAPMN[2], r1_GW_REVAP[2], r1_OV_N[2], r1_SLSUBBSN[2],
r2_CN2[2], r2_ESCO[2], r2_SOL_AWC[2],r1_ALPHA_BF[2], r2_GWQMN[2], r2_GW_DELAY[2],
r2_REVAPMN[2], r2_GW_REVAP[2], r2_OV_N[2], r2_SLSUBBSN[2],
r1_HRU_SLP[2], r1_USLE_K[2], r1_USLE_P[2], r1_P_UPDIS[2], r1_FRT_SURFACE[2],
r1_FRT_KG[2], r1_NPERCO[2], r1_ERORGN[2],
r2_HRU_SLP[2], r2_USLE_K[2], r2_USLE_P[2], r2_P_UPDIS[2], r2_FRT_SURFACE[2],
r2_FRT_KG[2], r2_NPERCO[2], r2_ERORGN[2],
r1_PPERCO[2], r1_CDN[2], r1_CMN[2], r1_PSP[2],
r2_PPERCO[2], r2_CDN[2], r2_CMN[2], r2_PSP[2])
par_init = c(-0.198073840010168,0.286267953999935,0.00746275760208839,
0.499157443870258,3500.6256700673,33.048081707854,339.59120477244,
0.032221701353556,0.514936793955879, -0.00241425477305834,0.0219494794735184,
0.41861314804868,-0.113031946972844, 0.69116116427558,1456.88879168975,39.7176103641318,
660.900428427947,0.0350548241647086, 0.559651736719884,0.380954113796417,
0.0156559405006919,0.0968484356612427,-0.0331742106779189,40.6194318342368,
0.386746499579996,-0.34870424019712,0.646162500043787,3.65240741158593,0.0113505467880799,
0.139763966103355,0.0545959533479474,50.2166504226717,0.337105246566447,-0.23731772825072,
0.886945397111459,3.03117861658487,13, 1.1,0.01,0.2, 13, 1.1,0.01,0.2)
# Optimización
head_log = c("FO","KGE23","KGE36", "KGE41",
"NSE23","NSE36", "NSE41",
"pbias23","pbias36", "pbias41",
"kge_PT_23", "kge_PT_41", "kge_NO3_41", "kge_NO3_23",
"pbias_PT_23", "pbias_PT_41", "pbias_NO3_23", "pbias_NO3_41",
"r1_CN2", "r1_ESCO", "r1_SOL_AWC","r1_ALPHA_BF","r1_GWQMN",
"r1_GW_DELAY","r1_REVAPMN","r1_GW_REVAP","r1_OV_N","r1_SLSUBBSN",
"r2_CN2", "r2_ESCO", "r2_SOL_AWC","r2_ALPHA_BF","r2_GWQMN",
"r2_GW_DELAY","r2_REVAPMN","r2_GW_REVAP","r2_OV_N","r2_SLSUBBSN",
"r1_HRU_SLP","r1_USLE_K","r1_USLE_P","r1_P_UPDIS","r1_FRT_SURFACE",
"r1_FRT_KG","r1_NPERCO","r1_ERORGN",
"r2_HRU_SLP","r2_USLE_K","r2_USLE_P","r2_P_UPDIS","r2_FRT_SURFACE",
"r2_FRT_KG","r2_NPERCO","r2_ERORGN",
"r1_PPERCO", "r1_CDN", "r1_CMN", "r1_PSP",
"r2_PPERCO", "r2_CDN", "r2_CMN", "r2_PSP")
fr = list.files("./run_optim_mensual/", pattern = "*.txt", full.names=TRUE)
file.remove(fr)
opt_pso = hydroPSO(par=par_init, fn = swat_optim,lower=par_lwr, upper=par_upr,
control = list(MinMax = "max",
maxit=500,
reltol=0.00001,
parallel="parallelWin",
par.nnodes=12,
par.pkgs = list("hydromad",
"SWATplusR",
"stringr",
"hydroGOF",
"tibble"),
write2disk=FALSE,
REPORT=1,
npart=48,
normalise=TRUE),
obs.data = obs.data,
project_dir = swat_project,
runrun_path = "./run_optim_mensual/",
head_log = head_log
)A work by Rafael Navas